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ABSTRACT 

We consider the dynamics of a protostellar disk surrounding a star in a circular-orbit binary system. 
Our aim is to determine whether, if the disk is initially tilted with respect to the plane of the binary 
orbit, the inclination of the system will increase or decrease with time. The problem is conveniently 
formulated in the binary frame in which the tidal potential of the companion star is static. We may then 
consider a steady, flat disk that is aligned with the binary plane and investigate its linear stability with 
respect to tilting or warping perturbations. The dynamics is controlled by the competing effects of the 
m — and m = 2 azimuthal Fourier components of the tidal potential. In the presence of dissipation, 
the to = component causes alignment of the system, while the to = 2 component has the opposite 
tendency. We find that disks that are sufficiently large, in particular those that extend to their tidal 
truncation radii, are generally stable and will therefore tend to alignment with the binary plane on a 
time-scale comparable to that found in previous studies. However, the effect of the to = 2 component 
is enhanced in the vicinity of resonances where the outer radius of the disk is such that the natural 
frequency of a global bending mode of the disk is equal to twice the binary orbital frequency. Under 
such circumstances, the disk can be unstable to tilting and acquire a warped shape, even in the absence 
of dissipation. The outer radius corresponding to the primary resonance is always smaller than the tidal 
| truncation radius. For disks smaller than the primary resonance, the to = 2 component may be able to 

cause a very slow growth of inclination through the effect of a near resonance that occurs close to the 
disk center. We discuss these results in the light of recent observations of protostellar disks in binary 
O 1 \ systems. 

Subject headings: accretion, accretion disks — binaries: close — hydrodynamics — instabilities — 
stars: pre-main sequence — waves 
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1. INTRODUCTION 

The existence of disks around young stars was spectacularly confirmed by direct images from the Hubble Space Telescope 
\ (HST) (McCaughrean & O'Dell 1996; Burrows et al. 1996). Observations suggest that young stars are usually found 
c/3 . in binary systems and that young binaries typically interact strongly with the disks that surround the stars (Ghez, 
■ Neugebauer, & Matthews 1993; Mathieu 1994; Osterloh & Beckwith 1995; Jensen, Mathieu, & Fuller 1996). There is 
growing evidence that disks within a binary are sometimes inclined with respect to the binary orbital plane. Such a case 
may have been seen in HST and Keck images of a disk in the young binary HK Tau (Stapelfeldt et al. 1998; Koresko 
1998). 

Suppose that a protostellar disk surrounds a star in a circular-orbit binary system, and that the disk is tilted with 
respect to the binary orbital plane. The evolution of the disk is affected by the tidal field of the companion star, as has 
been considered by Papaloizou & Terquem (1995). Some features of their analysis were confirmed in three-dimensional 
numerical simulations by Larwood et al. (1996). The basic physics involved may be summarized as follows (see also Bate 
et al. 2000). In a non-rotating frame of reference centered on the star about which the disk orbits, the companion star 
orbits at the binary frequency and exerts a time-dependent tidal torque on the disk. This torque may be decomposed 
into a steady component and an oscillatory component with a frequency of 2£!b, and their effects may be considered 
separately. 

Consider first the steady torque. If the disk were composed of non-interacting circular rings, the steady torque would 
cause each ring to precess, about an axis perpendicular to the binary plane, at a rate that depends on the radius of the ring, 
resulting in a rapid twisting of the disk. However, if the disk is able to maintain efficient radial communication, whether 
by wave propagation, viscosity, or self-gravitation, it may be able to resist this differential precession by establishing an 
internal torque in the disk. This can be arranged so that the net torque on each ring is such as to produce a single, uniform 
precession rate. However, to establish this internal torque, the disk must become warped. The concomitant dissipation 
changes the total angular momentum of the disk, tending to bring it into alignment with the binary plane in addition to 
causing accretion. 

Consider now the oscillatory torque. Applied to a single ring, this would cause a modulation of the precession rate 
and also a nutation (Katz et al. 1982). However, in the presence of radial communication, the oscillatory torque drives a 
bending wave (with azimuthal wavenumber to = 1) in the disk. Papaloizou & Terquem (1995) showed that, if the wave 
is subject to dissipation, it too may change the total angular momentum of the disk and tend to increase its inclination. 
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Tilting of protostellar disks 



The net effect of the steady and oscillatory torques determines whether an initially coplanar disk will acquire a tilt over 
time or whether an initially inclined disk will evolve towards coplanarity. The purpose of this paper is to determine this 
outcome, which could provide clues to the origin of misaligned disks in systems such as HK Tau. 

The basic mechanism suggested by Papaloizou & Terquem for generating a tilt by the oscillatory torque can be related 
to earlier work by Lubow (1992), who showed that an aligned, Keplerian disk in a circular binary may be linearly unstable 
to tilting if it contains a local resonance at which the orbital angular velocity £l(r) satisfies 



Here m* is the azimuthal wavenumber of the component of the tidal potential that is involved in the instability cycle. 
The cycle works through a mode-coupling process as follows: given a perturbation with m = 1 (a tilt), the tidal potential 
interacts with it to a drive a wave with to = to* — 1 at the resonant radius. This in turn interacts with the tidal potential 
to produce a stress with to = 1, which can influence the tilt. The role of dissipation is subtle, since some dissipation is 
required to provide a change in the angular momentum of the disk if instability is to occur, yet the associated damping 
can compete with the intrinsic growth rate of the instability. 

In particular, if the disk extends to the 3 : 1 resonance (Jl = 3f2b) it may be unstable to tilting through the to = 3 
component of the tidal potential. This resonance has the smallest to* for which equation (Q) can be satisfied (for a 
prograde disk) and is the closest resonance to the central star. It is difficult for disks to extend even as far as the 3 : 1 
resonance, because of the effects of tidal truncation (Paczyhski 1977; Papaloizou & Pringle 1977). Superhump binary 
disks might extend to the 3 : 1 resonance because of their extreme binary mass ratios, the secondary companion having 
less than 1/5 the mass of the primary about which the disk orbits (see the review by Osaki 1996). This instability is 
related to, and occurs at the same position as, the eccentric instability that is believed to be responsible for superhumps 
in cataclysmic variable disks (Lubow 1991). However, the growth rate is invariably much smaller for tilting than for 
eccentricity, and the weak tilt instability may be suppressed by the effects of viscous damping and accretion (Murray & 
Armitagc 1998). The same instabilities had been previously identified in the context of planetary rings for higher to* 
(Goldrcich & Trcmaine 1981; Borderies, Goldreich, & Tremaine 1984). More fundamentally, free-particle orbits undergo 
even stronger, parametric instabilities at these resonant locations (Paczyhski 1977), although free particles fail to model 
properly the behavior of a fluid disk at resonances. 

We relate this theory to the suggestion of Papaloizou & Terquem (1995) by noticing in equation (|l|) that, for to* = 2, a 
near resonance is obtained in the inner part of the disk where SI 3> fib- Indeed, Papaloizou & Terquem rely on the to = 2 
component of the tidal potential to drive an to = 1 bending wave in the tilted disk. The resulting response is a slowly 
rotating m = 1 bending wave, with frequency 2S7b in the inertial frame. Such a wave is close to resonance in the inner 
part of a nearly Keplerian disk because of the near coincidence of the effective wave driving frequency fl — 2f2b and the 
frequency of vertical oscillations Q z w f2; this is indeed the origin of equation ([j]) with m* = 2. An additional resonant 
effect occurs owing to the near coincidence of the driving frequency and the epicyclic frequency of horizontal oscillations 
K O. We describe the instability cycle associated with the oscillatory torque as a mode-coupling process in Fig. 1. 
However, because the resonance is not exact, and because of the importance of resonantly induced horizontal motions, a 
proper treatment requires a distinct analysis from that of Lubow (1992). 

In this paper, we therefore examine whether a flat, aligned disk in a binary is linearly unstable to tilting even if it 
does not extend to the 3 : 1 resonance. This problem is most conveniently analyzed in the binary frame where the tidal 
potential is static, since the disk can then be considered to be steady and to admit normal modes. These modes do not 
have a pure azimuthal wavenumber because the disk is non-axisymmetric as a result of tidal distortions. However, the 
tilting instability, if present, may be expected to appear as a modification of the rigid-tilt mode, which is trivial in the 
absence of the companion star. This mode may be followed continuously as the mass of the companion is increased, in 
order to determine whether it acquires a net rate of growth or decay. 

In general, the analysis of a normal mode of a tidally distorted disk is very difficult owing to the non-axisymmetric 
distortions of the disk. We therefore adopt the following simple approach, which is appropriate when only m = 1 bending 
waves are involved. We start by writing down the reduced equations for linear bending waves in a protostellar disk subject 
to an axisymmetric external potential (Section 2). These can be derived formally without great effort (see the Appendix). 
We then give a physical interpretation of these equations and use this insight to see how to modify them in the presence 
of a non-axisymmetric potential (Sections 3 and 4). We present a simple disk model (Section 5) and describe the results 
of numerical calculations of normal modes (Section 6). Some further analysis illuminates the underlying physics and helps 
to explain the numerical results (Section 7). Finally, we summarize our findings (Section 8). 



Consider a thin, non-self-gravitating disk in an external gravitational potential $(r, z) that is axisymmetric, but not 
necessarily spherically symmetric. Here (r, <f>, z) are cylindrical polar coordinates. The orbital angular velocity the 
epicyclic frequency /t(r), and the vertical frequency fl z (r) are defined by 3 



3 The true angular velocity of the fluid will depart from Q as a result of the radial pressure gradient and the vertical variation of the potential. 
Such departures generally depend on z and are of fractional order (H/r) 2 . They are fully taken into account in the analysis in the Appendix. 
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2. REDUCED DESCRIPTION OF LINEAR BENDING WAVES 



Lubow & Ogilvie 



3 



K 2 = 4tt 2 + 2rf^, 
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We consider a situation in which the disk is nearly Keplerian and almost inviscid in the sense that 
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where H(r) is the semi-thickness of the disk and a the dimensionless viscosity parameter. Then the linearized equations 
for bending waves (with azimuthal wavenumber m = 1) may be written 
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Here S(r) is the surface density and I{r) the second vertical moment of the density, defined by 

£= / pdz, 1=1 pz 2 dz. 
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The dimensionless complex variable W(r, t) is defined by W = £ x + i£ y , where £(r,t) is the tilt vector, a unit vector 
parallel to the local angular momentum vector of the disk. The complex variable G(r, t) represents the internal torque 
which acts to communicate stresses radially through the disk (see below). 

The derivation of these equations may be found in the Appendix. Equivalent equations, although presented in quite 
different notations, have been derived by Papaloizou & Lin (1995) and Demianski & Ivanov (1997). In addition to 
conditions (|))-(0), it is required that the warp vary on a length-scale long compared to the thickness of the disk, and on a 
time-scale long compared to the local orbital time-scale. However, any evolution of the disk on the (much longer) viscous 
time-scale is neglected. The (dynamic) viscosity is assumed to be isotropic and proportional to the pressure (p = ap/il). 

We emphasize the physical interpretation of these equations. Equation (||) contains the horizontal components of the 
angular momentum equation encoded in the combination l x + iy'. In vectorial form it may be written 



The second moment is related to the integrated pressure through 

pdz = 2X2 2 , . 



Zr 2 n 



di 
dt 



r dr 



(12) 



where 2TrG(r,t) is the internal torque and T(r,t) the external torque density acting on the disk. In the present case the 
external torque arises from a lack of spherical symmetry in the potential. The complex variable G is simply G x + iG y . 
Its equation may be written, in vectorial form, 
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The internal torque is mediated by horizontal epicyclic motions that are driven near resonance by horizontal pressure 
gradients in the warped disk, an effect identified by Papaloizou & Pringle (1983). The horizontal motions are proportional 
to z and are therefore subject to strong viscous dissipation which is the dominant channel of damping of the bending 
waves. 

We note that slowly varying m = 1 bending waves or warps may be quite generally described by conservation equations 
for mass and angular momentum (Pringle 1992; Ogilvie 1999). The relevant relation of G to I and its derivatives, however, 
depends strongly on the thickness of the disk, the viscosity, and the rotation law. In this paper we are considering a 
parameter range appropriate to protostellar disks, and are assuming that the warping is small so that a linear theory is 
valid. For disks in which a £ (H/r), see Papaloizou & Pringle (1983), Pringle (1992), and Ogilvie (1999, 2000). 
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3. TIDAL TORQUE ON A TILTED RING 

There are two dynamical degrees of freedom in the system described by the above equations. One is the tilting of 
the disk at each radius according to the tilt vector i. The other is the horizontal motions described by G, which cause 
eccentric distortions of the disk that are proportional to z. In spite of this complexity, the external torque density T in 
equation (|l2|) is very simple (cf. eq. || ) : it is equal to the torque exerted by the external potential on a disk composed of 
arbitrarily thin circular rings of uniform density that are tilted according the tilt vector £. The eccentric distortions may 
be disregarded when calculating the external torque to the required order. ^ 

We proceed to derive an expression for the torque exerted by the full potential of the companion star (of mass M 2 ) 
on a tilted ring of the disk, treated as a thin and narrow circular ring of radius r and uniform density. Adopt Cartesian 
coordinates (x, y, z) with origin at the center of the ring, and with the ring in the xy-plane. Then the position of an 
arbitrary point on the ring is 

r = (r cos 4>, r smcj), 0), (14) 

where <f> is the azimuthal angle measured around the ring. Assume, without loss of generality, that the companion star 
lies instantaneously in the ccz-plane at position 

r h = (r b cos/3, 0,r b sin/3), (15) 
where r b is the binary radius and f3 the angle of inclination. Then the force per unit mass at position r on the ring is 

GM 2 (r h -r) 
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and the corresponding torque per unit mass is 
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The azimuthally averaged torque per unit mass is 
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The x- and z-components vanish owing to the antisymmetry of the integrands. The remaining component is 
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In general, this may be expressed in terms of elliptic integrals. For small (3, however, we have 
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where 



is the Laplace coefficient. In vectorial form, therefore, the torque density is 
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with fractional corrections of 0(W 2 ). 



4. DYNAMICS IN THE BINARY FRAME 



We now consider the dynamics in the binary frame, which rotates with angular velocity Jl b = fi b e z . This requires that 
we replace, in equations ( |T2"| ) and (|T3|), 
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With the companion star located on the positive a;-axis, we have, in linear theory, 
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4 The effects of eccentric and tidal distortions are considered implicitly in Section 7, where they are found to be unimportant for the linear 
growth rates we derive. 
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The orbital, epicyclic, and vertical frequencies are all calculated using the m = total potential. This gives 
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where M\ is the mass of the star about which the disk orbits. We assume, without loss of generality, that fl > 0, but 
allow for the orbit of the companion star to be either prograde or retrograde according to 



The final equations are 



£r 2 ft 



± 



dL 



G{M X + M 2 ) 



1/2 



at 



l8G r 



r dr 



Sr 2 n (M± + n h i x 

dG x GM 2 
dG 



1 dG y GM 2 
r dr 



at 



+ QbG x 



ir^rQ 
GM 2 



J 3/2 



b?l [ - 



Ha 



3/2 

G v + anG. 



J 3/2 



n, 



G, 



*fiG„ = 



Er4, 

lr 3 n 3 d£ x 
4 ~dr' 

lr 3 fl 3 d£ y 
4 ~dr' 



(28) 

(29) 
(30) 
(31) 
(32) 



Since the coefficients of these equations are independent of time, we may seek normal modes of the form 

£ x (r,t)=Re[t x (r)e iut ] , (33) 

£j / (r,t)=Re[4(r)e iw *] , (34) 

G x (r,t)=Re[G x (r)e wt ] , (35) 

G y (r,t) =Re [G,(r)eM , (36) 

where u> is a complex frequency eigenvalue. The problem has then been reduced to solving an eigenvalue problem involving 
a fourth-order system of ordinary differential equations (ODEs). 

If we return to the original complex notation, we find that the equations have become non-analytic, effectively increasing 
the order of the dynamical system: 
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iG + aQG = 
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In the combination W + W* , the term W arises from the m — component of the tidal potential (cf. eq. |§), while the 
non- analytic term W* arises from the m = 2 component. In the normal- mode solution, W and G have the form 
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The motion seen in the inertial frame is more complicated than a single mode. We have 

£. £ x & x -\- £y Gy -\- £% (^5) 

where (e x , e y , e z ) are unit vectors in the binary frame. These are related to the unit vectors (e x ,e y ,e z ) in the inertial 
frame by 

e x - ie y = (e x - ie y ) e i0bt , (46) 



and so 
where 



^x &x ~t~ ^y &y ^-z ^z: (^7) 



(48) 



Here u> = w r + zcjj. Therefore two components are seen in the inertial frame, which have distinct frequencies, \u) T ± £!b|, 
but the same rate of growth or decay. 

5. DISK MODEL 

For simplicity, we assume that the vertical structure of the unperturbed disk is that of a polytrope of index n. To 
satisfy vertical hydrostatic equilibrium, the density distribution, for a thin disk, is then of the form 

(2 \ n 

where H(r) is the semi-thickness. The surface density and second moment are related by 

1=7^- (50) 
2n + 3 v J 

For the radial structure, we specify 

H , s 

- = e (51) 
r 

and 

E = Eor- 1 / 2 /, (52) 

where e is a small constant, Eo an arbitrary constant, and f(r) a function that is approximately equal to unity except 
near the inner and outer radii of the disk, where it tapers linearly to zero. Over most of the disk this gives approximately 
H oc r, E oc r -1 / 2 , and X oc r z / 2 . 
For the tapering function, we take 

/ = tanh ( — — ) tanh ( , (53) 

\ Wi J \ w 2 J 

where T\ and r 2 are the inner and outer radii of the disk, and W\ and w 2 are the widths of the tapers near each edge, which 
are taken to be equal to the local semi-thickness. With / tapering linearly to zero, the edges are regular singular points 
of the governing equations. The appropriate boundary condition in each case is that W should be regular there, which 
implies that G vanishes. Clearly the internal torque cannot be transmitted across a free boundary of the disk. However, 
if the inner disk were terminated by a magnetosphere, for example, this boundary condition may require modification. 

This model is very similar to that used by Papaloizou & Terquem (1995) except that the disk has an inner edge. For 
reasons that we explain in Section 7, we do not attempt to impose an 'ingoing wave' boundary condition at the center of 
the disk. 



6. NUMERICAL RESULTS 

Equations (|2^)-(^) are solved numerically using the complex variables defined in equations (|3^)-(|3^). When solving 
the ODEs for a normal mode, it is advisable to integrate away from the singular points at the edges of the disk. We 
apply the arbitrary normalization condition £ x (^i) = 1 an d guess the values of the four complex parameters u), £ y (ri), 
Vxfo), and t y (r2)- We then integrate separately into r > n and r < r%, meeting at the midpoint, where £ x , £ y , G x , and 
G y should all be continuous. These four conditions are solved by Newton-Raphson iteration, using derivative information 
obtained by simultaneously integrating the ODEs differentiated with respect to the four parameters. 
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Table 1 

Parameters of the reference model. 



Parameter 


Symbol 


Value 


Mass ratio 


q = M 2 /Mi 


1 


Angular semi-thickness 


€ 


0.1 


Viscosity parameter 


a 


0.01 


Inner radius 


n/r h 


0.01 


Outer radius 




0.3 


Width of inner taper 


wi/n 


0.1 


Width of outer taper 


w 2 /r 2 


0.1 


Polytropic index 


n 


3/2 



Table 2 

Frequency eigenvalues for the reference model, but with q = and a = 0. 



Mode 


ti>/Ob 


Mode 


w/fib 





-1 


0* 


1 


1+ 


-0.3074 


i; 


+0.3074 


1_ 


-1.6926 


1* 


+ 1.6926 


2+ 


+0.3163 




-0.3163 


2_ 


-2.3163 


2*_ 


+2.3163 


3+ 


+0.9213 


y + 


-0.9213 


3_ 


-2.9213 


3* 


+2.9213 


4+ 


+ 1.5180 


A* + 


-1.5180 


4_ 


-3.5180 


4* 


+3.5180 



Table 3 

Frequency eigenvalues for the reference model, but with q = and a = 0.01. 



Mode w/Ob 






-1 


1+ 


-0.3101 + 0.0752i 


1_ 


-1.6899 + 0.0752i 


2+ 


+0.3144+ 0.0925* 


2_ 


-2.3144 +0.0925i 


3+ 


+0.9197 + 0.1020i 


3_ 


-2.9197 + 0.1020i 


4+ 


+1.5166 + 0.1087z 


4_ 


-3.5166 + 0.1087z 



8 



Tilting of protostellar disks 



6.1. Reference model 

We first identify a 'reference model' with parameters that we consider appropriate for a protostellar disk that is tidally 
truncated by the companion star (Table 1). The orbit of the companion is taken to be prograde. 

Before considering the reference model as such, we examine the same disk but with no viscosity (a — 0) and with a 
companion of zero mass (q = 0). An infinite set of discrete bending modes is obtained, which are characterized by the 
number of nodes in the eigcnfunction l x (say). The basic frequencies of these modes in the inertial frame are loq = 0, 
ui\ = 0.6926 0b, ^2 = 1.3163 fib, &3 = 1.9213 fib, uj 4 = 2.5180 fib, etc. We refer to these modes as modes 0, 1, 2, 3, 4, etc. 
Mode is the (trivial) rigid-tilt mode and has no nodes. 

In the binary frame, the full set of frequencies appears much more complicated, as shown in Table 2. The modes in 
the left-hand column consist purely of W+ and G+, having W- = and G_ = 0. For such a mode, the frequency in the 
binary frame is less than the frequency in the inertial frame by an amount fib (cf. eq. p8[). Since, in the inertial frame, 
we may have a prograde or retrograde mode n with frequency ±w„, we obtain frequencies ico n — fib in the binary frame. 
These are labeled n±. The modes in the right-hand column are physically equivalent. The eigenfunctions and eigenvalues 
are obtained from those in the left-hand column by complex conjugation and a change of sign. Such modes consist purely 
of W- and G_, having W+ = and G+ = 0. Thus the frequencies in the binary frame are =Fcj„ + fib- These modes are 
labeled ri±. 

We consider next the effect of a small viscosity on the modes by increasing a from to its reference value 0.01, but still 
with a companion of zero mass (q = 0). The results are shown in Table 3. We omit the complex-conjugate modes from 
now on, but their existence should not be forgotten. Evidently the real part of the frequency changes very little in the 
presence of a small viscosity, but, with the exception of the rigid-tilt mode, the frequency acquires a positive imaginary 
part, which signifies a damping rate. The damping rate depends relatively little on the order of the mode. It can be seen 
from the governing equations that the effect of viscosity is simply to damp the horizontal motions locally at a rate afi 
(cf. eq. Jig]). Since the horizontal motions are an essential part of each proper bending mode, this leads to a damping 
rate for each mode of order afi (evaluated in the outer parts of the disk). The exception is the rigid-tilt mode, for which 
the horizontal motions are exactly zero. 

Finally, we reach the reference model by increasing the binary mass ratio q from to its reference value 1. We start 
with mode 0, which corresponds to a rigid tilt and consists purely of W+. The frequency of the mode (now the 'modified' 
rigid-tilt mode) changes continuously from —fib to (—1.0484 + 0.000258i)fib- The mode also acquires a W_ component. 
Viewed in the inertial frame, the mode changes from a pure W+ mode with zero frequency to a combination of W+ and W- 
contributions having frequencies of 0.0484 fib and 2.0484 fib, respectively (see eq. Q). The first frequency corresponds to 
a retrograde precession of the tilted disk, forced by the m = component of the tidal potential. The second corresponds 
to the forcing of a bending wave (W_) by the m = 2 component of the potential. The two potential components provide 
the 'steady' and 'oscillatory' torques, respectively. Since the imaginary part of the frequency is positive, the whole pattern 
decays at a rate 0.000258 fib- The other modes of the disk are of course damped much more rapidly, and we conclude 
that the reference model disk is linearly stable to tilting. 

6.2. Resonances 

We now search the parameter space around the reference model for any regions of instability. In particular, we try 
varying the outer radius of the disk. In Fig. 2 we plot the dimensionless growth rate — Wj/fib against r^/r^ for a number 
of different values of a. It is clear that the net growth rate is a combination of two parts. One part is a damping (uj\ > 0) 
that is proportional to a and increases rapidly with increasing T2 ■ The second part is a growth (wj < 0) with an entirely 
different behavior. The growth is localized in a sequence of peaks which become higher and narrower as a decreases. In 
Fig. 3 we show an expanded view of the primary peak for the cases a = and a = 0.001. 

To verify the origin of the two parts, we repeated the calculation using equations that retain only the m = component 
of the tidal potential, or only the m = 2 component. It is obvious from this that the damping is due entirely to the m = 
component, while the growth is due entirely to the m = 2 component. There is a slight shift in the positions of the peaks 
when the m = component of the tidal potential is neglected. 

It is evident that the growth (that is, the instability) is associated with a series of resonances that occur when the outer 
radius of the disk is in the vicinity of certain discrete values. In the absence of viscosity, the resonances come about as 
follows. As T2/rb is varied, the frequency eigenvalues of all bending modes migrate along the real axis in the cj-plane. 
With the exception of mode 0, all modes are very sensitive to the position of the outer boundary, which reflects the 
waves. As a result, collisions occur on the real axis. In particular, when ^/r^ is increased from 0.1 towards the primary 
resonance, mode undergoes a collision with mode 11 (a bending mode with one node). The modes move briefly off 
the real axis, producing a complex-conjugate pair, and then return to the real axis to continue their original migration. 
The other resonances occur when mode undergoes collisions with modes 23_, 35_, etc. During a collision, the two modes 
exchange characteristics, and the eigenfunctions are hybrids of the two original ones. In particular, mode no longer 
resembles a rigid tilt during a collision with a proper bending mode. This means that a disk made unstable by this means 
would develop a warped shape (see Section 6.4 below). 

In the presence of a very small viscosity, the proper bending modes are damped and their eigenvalues are displaced 
somewhat above the real axis. The collisions are no longer exact and each mode can be followed continuously as ^/ri, is 
varied. For a = 0.001, say, the modes pass sufficiently close that a strong interaction occurs. The tracks of the eigenvalues 
are deflected to avoid a collision, and, in so doing, mode acquires a positive growth rate that appears as a resonance. 
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During the interaction, the eigenfunction of mode is distorted significantly from a rigid tilt, but not so strongly as in 
the inviscid case (see Section 6.4 below). 

When the viscosity is increased, the resonances become broader and weaker. A positive growth rate is not achieved 
if the height of the resonance is less than the damping rate arising from the m = potential. Therefore the regions of 
instability are suppressed as a is increased. It appears that, as long as the primary resonance survives, the net growth 
rate (for a > 0) is always positive for disks smaller than the size of the primary resonance, although the growth rate may 
be minuscule. This may be considered as a long tail of the primary resonance. However, the primary peak is dramatically 
reduced in height as a is increased, and it also shifts to smaller radius. In the cases investigated here, all traces of 
instability are eliminated when a = 0.1. 

To elucidate further the condition for resonance, we examined the bending modes at their points of collision with mode 
and evaluated their natural frequencies (i.e. in the absence of the tidal potential, and evaluated in the inertial frame). 
In each case the natural frequency is close to 2f2b at the point of collision. The resonances therefore occur when the 
oscillatory torque due to the m = 2 potential resonates with a free bending mode of the disk. 

We remark that the global resonant excitation of bending waves has been identified by Terquem (1998) when calculating 
the tidal torque exerted on a protostellar disk by a companion in an inclined circular orbit. However, the consequences 
for the evolution of the relative inclination of the system were not investigated. 

The results for a companion in a retrograde orbit are not significantly different. The heights of the resonant peaks are 
very similar, but they are shifted slightly in radius. The shift of the resonances (also observed, as noted above, when the 
m = component of the tidal potential is omitted) is related to the precession of the disk, which changes the effective 
frequency of the oscillatory torque and, therefore, the condition for resonance. The precession is always retrograde in the 
inertial frame, irrespective of the sense of the companion's orbit. Therefore the effective driving frequency depends on the 
sense of the orbit, but the shift is generally small. 

6.3. Precession rate and decay rate 

In Fig. 4 we plot the precession rate of the modified rigid-tilt mode against the outer radius of the disk, for the reference 
model. The precession is always retrograde and the rate increases rapidly with increasing r-i. Excellent agreement is 
found with the simple analytic approximation given by Bate et al. (2000; eq. [22]). (We have set the dimensionless 
parameter K — 0.4, since this represents fairly accurately the disk models we have adopted.) For much smaller values 
of a, a noticeable deviation from this curve occurs in the vicinity of resonances, since the path of the eigenvalue in the 
w-plane is temporarily diverted. 

In Fig. 5 we plot the decay rate of the modified rigid-tilt mode, for the reference model. When only the m = component 
of the tidal potential is included, the decay rate is always positive and increases rapidly with increasing r^- When the full 
potential is used, the behavior is modified in the vicinity of resonances. Also shown is the simple estimate 1/ialign given 
by Bate et al. (2000; eq. [35]). Apart from the resonances, the simple estimate captures the correct dependence on r-i. It 
should be borne in mind that the estimate of Bate et al. (2000) was based on an order-of-magnitude analysis, and can be 
expected to be accurate only within a factor of order unity. 

6.4. Shape of the disk 

For comparison with observations, it is of interest to examine the shape adopted by the disk while executing the 
modified rigid-tilt mode. Information on the shape of the disk is contained in four real functions of radius, namely the 
real and imaginary parts of the eigenfunctions £ x (r) and £ y (r). We display this information in Figs 6 and 7 by showing 
cross-sections through the disk in the xz- and yz-planes at two instants, corresponding to phase and phase 7r/2 of the 
period seen in the binary frame. 

Fig. 6 is for a disk with r^jr^ = 0.118, in the middle of the primary resonance. Three different viscosities, a = 0, 
a = 0.001, and a — 0.01, are considered. In each case the mode has a positive growth rate. In the absence of viscosity, 
the resonance is strong and the disk becomes distinctly warped in a smooth and global manner. As already noted, when 
viscosity is included, the resonance is much weaker and the disk appears tilted with less noticeable warping. 

Fig. 7 is for a disk with the reference value r 2 /rb = 0.3 representative of a tidally truncated disk. We fix a — 0.01 and 
consider disks of varying thickness, e = 0.1, e = 0.05, and e = 0.03. In each case the mode is damped. For e = 0.1, the 
disk appears tilted without noticeable warping. For thinner disks, the deviation from a rigid tilt is noticeable in the outer 
part of the disk where the tidal forcing is strongest. 

Recall that the derivation of equations ( |l2| ) and (|l^) requires that the warp vary on a length-scale long compared to 
the thickness of the disk (see the Appendix] - ^ This condition is indeed satisfied in the solutions we present here. 

7. EXPANSION IN THE TIDAL POTENTIAL 
7.1. Basic equations 

The normal-mode description affords an especially compact representation of the dynamics and is very suitable for the 
numerical analysis. In this section we 'unpack' the eigenfunction to reveal the essential physics of the problem. We write 
the basic equations in the general form 

Zr 2 n(^+in h w) = -^--i(AW + BW*), (54) 
\ at J r or 
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— — I- iflhG — i(CG + DG ) + atiG — — — , (55) 

ot 4 or 

with unspecified coefficients A, B, G, and D arising from the tidal potential. In view of our earlier discussion, terms A 
and C are due to the m = component of the potential, while the non-analytic terms B and D are due to the m = 2 
component. We allow for the possibility that tidal distortions of the disk may introduce additional complexities (such as 
a term D) that we have not foreseen. 

For a normal mode of the form (]3^)-(f40|), we have 

1 dG 

(iuu + iVL h )Y>r 2 VLW+ = -± - i(AW+ + BW1), (56) 

r dr 

{-iu* + if2 b )£r 2 fiW_ = - i(AW- + BW1), (57) 

r dr 



Tr 3 fi 3 dW, 

(iuj + in h )G+ - i(CG+ + DG*_) + a nG+ = -±, (58) 

4 dr 

2> 3 3 dW 

i-iu>* + in h )G_ - i(GG_ + DG*,) + aftG_ = — . (59) 

T A dr 

7.2. Expansions 

We now expand the equations in powers of the tidal potential, indicated by a numerical superscript. The unspecified 
coefficients may be assumed to have expansions 

A = A^ + A^ + • • ■ , (60) 
etc., since they vanish in the absence of the tidal potential. The eigenvalue and eigenfunction have the expansions 





u w + ..., 


(61) 


w+ = + 




+ wf> + ..., 


(62) 


W- = 


w w 


+ w^ + ---, 


(63) 


G+ = 




+-0?+-, 


(64) 


G_ = 


G w 


+ G (2) + ---, 


(65) 


rigid-tilt mode with 








^ (0) = -n b , 


w 


1°^ = constant. 


(66) 



The rigid-tilt amplitude could be arbitrarily specified as wi° = 1 , but we retain for clarity in the equations below. 



At first order, we obtain 



7.3. Solution 
-(i) 



r dr 



(67) 



ldG W 



2in h Er 2 QWP = - l B^wf ) \ (68) 

r dr 

aQG W = _ (69) 

m 2r 3 n 3 dW (1) 

(2in h + att)G W = ~ . (70) 

4 dr 

From equation (|67j), using the fact that G vanishes at the edges of the disk, we immediately obtain the solvability condition 

cjW = - A^W^rdr/ & 2 fW{ 0) rcZr, (71) 

Jr\ I Jri 

which relates the precession rate (at first order) to the total horizontal tidal torque on the disk (at first order) divided by 
the horizontal angular momentum of the tilted disk. This effect is due to A and therefore to the m = component of the 
tidal potential. Equations (67)-j7C|) can then, in principle, be solved for W± and G£ . 
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At second order, we obtain 

r dr \ / 



plus three further equations, which will not be required. This time the solvability condition is 

f»"2 



to 



.2) 



^ U^wf + > + u;«Er 2 QWf 5 + B^W W *) rdrl C 

i I Jri 



T,r 2 QW^ ) rdr. 



If, as we assume, A is real, then a/ 1 ) is real. After some further manipulations we then obtain 



Im 



( Aa 



\lr 4 tt 



G 



(i) 



G 



(i) 



r dr 



(0) 



r dr. 



(72) 



(73) 



(74) 



This shows that G} , which is caused by the m = component of the potential, causes pure damping, while G_ , which 
is caused by the m = 2 component, causes pure growth. The net effect depends on which is larger in the norm defined 
above. Note that coefficients C and D have no effect to this order. Also, the second-order coefficient A^ 2 \ which we did 
not attempt to calculate, does not affect the growth or decay rate at second order (although it does affect the precession 
frequency at second order). Such a coefficient could arise because the tidal torque on a tilted ring may have a second-order 
correction owing to the tidal distortion of the ring. 

Furthermore, since is independent of a according to equation (|67]), the damping is simply proportional to a. The 

dependence of the growth on a is less clear since G^} itself depends on a in a complicated way according to the coupled 



equations (pq) and (|7C|). Some insight into these equations is obtained by considering the case a — 0, for which we find 



d 2 W { }> d\n{lr 3 n 3 ) dW_ 
dr 2 



r(l) 



dr 



dr 



lggng (i) 8Q h BWwj oh 

m 2 



ir 2 n 3 



(75) 



For our disk model, we have approximately f2 oc r 3 / 2 , Ear 1 / 2 , and I cx r 3 / 2 over most of the disk. We then obtain, 
approximately, 



d 2 W_ 



(i) 



where 



dr 2 



X 



-W (1) = 

A 3 VV - 



16(2n + 3)(l + g) 



ir 2 n 3 

-,1/3 



(0)* 



(76) 



(77) 



and e is the angular semi-thickness H/r of the disk. This is an inhomogeneous Airy equation such as is common in 
problems of resonant wave excitation in differentially rotating disks. Here, the resonance is at the exact center of the disk, 
in accord with equation ([!]). The forcing term on the right-hand side, however, is proportional to r 5 / 2 and is therefore 
concentrated in the outer parts of the disk. 

7.4. Interpretation 

The magnitude of the response (and therefore G^) depends on the overlap between the forcing function and 

the solutions of the homogeneous equation, Ai(—r/X) and Bi(—r/X). One may then distinguish two cases, depending on 
whether the outer radius r2 satisfies r2 > A or r2 ^ A. 

If r-2 ^> A, the homogeneous solutions are highly oscillatory over the disk and the overlap will be very small unless a 
global resonance occurs. This happens when there is a homogeneous solution that (nearly) satisfies both radial boundary 
conditions. This means, in fact, that the frequency of a free bending mode of the disk (in the inertial frame) is (nearly) 
equal to 2f2b- Then the operator on the left-hand side of equation (f7q) is (nearly) singular and a large response results. 
This clearly occurs during the inviscid resonances. When this happens, equation (J7q) breaks down; however, the analysis 
in Section 6 is valid. 

If, instead, r2 <^ A, the response is also reasonably large because there is little or no cancellation in the overlap integral. 
This can explain the long tail of the primary resonance, where the net growth rate is found to be positive for sufficiently 
small, but non-zero, a. 

For the reference model, A « 0.037rb- As r 2 /rb is reduced from 0.3 to 0.05, we pass from the first case r 2 ^> A, through 
several resonances, towards the second case, r-i ~ A. The interpretation given above can therefore explain the behavior 
found in Section 6. 

It is natural to ask whether the parameters of real disks are likely to allow a tilting instability in practice. We have 



tested how the the value of the outer radius at which the primary resonance occurs, r 2 



varies with all the parameters 



of the model. The variations with e, n, and q are well approximated by 



3.2A, A being given by equation U77J). The 
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variations with ri/rb, ioi/ri, and W2A2 are all less significant. Comparing this estimate of r p with the tidal radius r t of 
a disk estimated by Papaloizou & Pringle (1977), we find that the inequality 

^"(sf 3 ' < 78 » 

is satisfied for 0.2 < q < 10. We conclude that tidally truncated disks extend too far beyond the primary resonance for 
instability to occur, unless e ;> 0.4, which is not suggested by observations (although it should be remembered that H 
is the true semi-thickness of our polytropic model, and not an approximate scale-height). For tidally truncated disks, a 
tilting instability would occur only in the unlikely case that a higher-order resonance condition were met. 

It can also be seen from the above that to impose an 'ingoing wave' boundary condition at the center of the disk, as was 

done by Papaloizou & Terquem (1995), is questionable. Those authors envisaged that the bending wave (i.e. W^) would 
be excited at the outer edge of the disk and would propagate inwards, growing in amplitude until nonlinear effects caused 
it to dissipate. The wave would then fail to reflect from the center of the disk. In contrast, we find that the wave may 
be considered to be launched at a local resonance located at the center of the disk. Since the tidal forcing vanishes there, 
the wave is not significantly excited unless the width of the resonance (proportional to A) becomes comparable to the 
radius of the disk. In this case, the wave is launched at all radii and global resonant effects must be taken into account. 
However, the wave amplitude does not diverge at the center of the disk; equation (|7^) has a solution with a finite tilt and 
vanishing torque at r = for the surface density profile adopted. Nonlinear dissipation may not occur. Indeed, perhaps 
contrary to conventional wisdom, the instability at the primary resonance (where the disk edge is at radius r p ) operates 
in a completely inviscid disk with reflecting boundaries. 

The effects of the contribution of W_ to the tilt growth can be related back to the mode-coupling description seen 
in Fig. 1. In particular, wave equation (|7^) describes the generation of the wave W- through the driving term on the 
right-hand side of that equation. This term involves the interaction of the m = 2 tidal potential, represented by with 

the rigid tilt, W+ . This interaction produces a wave, W- , of the form of an m — 1 bending wave having frequency nearly 
equal to 217b in the inertial frame. The interaction of the wave with the tidal field produces a stress that corresponds to 
the tilt growth-rate contribution B^W_* in equation (73). Therefore the instability mechanism described by Lubow 



(1992) is always at work here, but with the differences that global resonant effects can be important, and that dissipation 
is not required. Furthermore, our comparison with the damping rate induced by the m — potential indicates that the 
instability is suppressed for tidally truncated disks, except in the unlikely event of a high-order resonance. 

8. SUMMARY AND DISCUSSION 

In this paper, we have considered the linear stability of a coplanar protostellar disk that surrounds a star in a circular- 
orbit binary system. We have determined whether a slight tilt introduced into the disk would grow or decay in time. The 
outcome depends on the size of the disk. For disks that are truncated by standard tidal torques, typically resulting in an 
outer disk radius of about 0.3 times the binary separation, we find that the disk tilt generally decays in time. 

For smaller disks, tilt growth is possible. As seen in Fig. 2, a disk undergoes a strong, 'primary' resonance with the 
tidal field when its outer radius is a certain fraction of the binary separation. This characteristic radius, which we denote 
by r p , is approximately 0.118 times the binary separation for the parameters we have considered (see Table 1), but would 
be smaller still for thinner disks with H/r < 0.1. In such a resonance, the disk experiences a growing tilt and becomes 
significantly warped (see Fig. 6). This resonance occurs when the frequency of the lowest-order global bending mode in 
the disk matches the tidal forcing frequency, which is here twice the binary orbital frequency. Weaker resonances occur 
at a series of discrete resonances corresponding to radii greater than r p . There is also a near resonance that occurs close 
to the disk center. For disks smaller than radius r p , this resonance causes a very slight tilt growth if a is sufficiently small 
(but non-zero), and any initial tilt would be retained. 

For disks with radii larger than r p , including disks truncated by standard tidal torques, the tilt will decay on approx- 
imately the viscous time-scale of the disk, or roughly 10 3 binary orbits for a = 0.01 (see Fig. 5). For disks with large 
tilts, nonlinear effects may shorten the time-scale to reach small tilts, perhaps to the precessional time-scale of the disk, 
or about 20 binary orbits (Bate et al. 2000). 

The net outcome of growth or decay of the disk tilt is determined by the competition of two torques. As seen in the 
inertial frame, the tidal torque acting on a tilted disk may be decomposed into a steady component and an oscillatory 
component with twice the binary orbital frequency. The steady torque, resulting from the m = component of the tidal 
field, causes the disk to become aligned with the binary orbit in the presence of dissipation, while the oscillatory torque, 
resulting from the m = 2 component of the tidal field, causes misalignment. The steady torque produces an intuitively 
simple result because it causes the disk to settle to a state of coplanarity, where it experiences a minimum tidal potential 
energy, as a result of dissipation. The effect of the oscillatory torque is somewhat counterintuitive, but can be understood 
in terms of a mode-coupling model (see Fig. 1). Provided that a is sufficiently small, the oscillatory torque slightly 
dominates for smaller disks because material in such disks is generally closer to the near resonance that occurs in the 
vicinity of the disk center (see eq. [|l|). 

A major issue is the origin of the tilt in observed protostellar disks. In the case of HK Tau, the disk surrounds the 
secondary star, but the two stars are similar in spectral type (Monin, Menard, & Duchene 1998). Although there are 
considerable uncertainties in the system parameters, the disk could extend to its standard tidal truncation radius, as 
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suggested by Stapelfeldt et al. (1998). In that case, the results of this paper imply that tidal effects may cause decay of 
the primordial tilt, but in any case would not cause tilt growth. On the other hand, the existence of the tilt means that 
the decay time-scale cannot be much shorter than the binary age, estimated as 5 x 10 5 yr. This places some constraints 
on both the theory and the binary parameters, although there are considerable uncertainties. For example, consider the 
case that the binary separation is close to its projected value of 340 AU. For a — 0.01, the linear tilt decay time-scale 
(based on Figure 2) would be several times longer than the estimated system age. On the other hand, the nonlinear decay 
time-scale estimate of Bate et al (2000) suggests a decay time-scale substantially shorter than the estimated age. The 
nonlinear time-scale estimate would be more compatible with a somewhat larger binary separation. 

The predicted shape of a tilted, tidally truncated disk with H/r w 0.1 is not strongly warped (see Fig. 7), in accord 
with the observations (Stapelfeldt et al. 1998; Koresko 1998). The lack of an observed warp cannot be used as evidence 
against binarity. On the other hand, a slight warp does occur for thinner disks such as in Fig. 7 case b, which could be 
observed as a small asymmetry. 

Note that the decay time-scale of proper bending modes of this disk (based on Table 3) is of order 10 4 yr if a = 0.01, 
much shorter than the linear tilt decay time-scale. If the disk were tilted and warped in an arbitrary way as a result of 
its formation process, we would expect it to evolve rapidly to a tilted but essentially unwarped shape, then the tilt itself 
would decay on a longer time-scale. However, the nonlinear effects discussed by Bate et al. (2000) are likely to speed up 
both stages considerably. 

Similar considerations apply to a recent observational test of coplanarity among a sample of T Tauri binaries by Donar, 
Jensen, & Mathieu (2000). The data show some evidence for approximate coplanarity on a statistical basis. It is possible 
that some tidal evolution of the tilt towards coplanarity may be have occurred, if the tilt decays as rapidly as a disk 
precessional time-scale. 

It is important to understand whether disk truncation could occur close to the resonant radius r p , so that the disk 
would be unstable to tilting. For disks with a substantial tilt, Terquem (1998) has shown that a disk of radius close to 
r p is sometimes subject to a strong resonant torque that is parallel to its spin axis. This resonant torque can exceed the 
viscous torque in the disk for sufficiently small values of a, a <; 10~ 3 . If this torque could truncate an initially tilted disk 
at radius r p , the disk might become strongly warped (as seen in Fig. 6) and tilted further. The disk radius would be less 
than half of the standard tidal truncation radius. However, it is unclear that this torque would lead to disk truncation at 
r p because it is smoothly distributed over the disk rather than being concentrated near r p . This is because the resonance 
is global rather than local. The lack of a strong warp in HK Tau argues against this process in that system. 

Disks in cataclysmic binaries are expected to be much colder than protostellar disks, having a smaller value of H/r. 
Consequently, such disks are even less likely to be unstable to tilting as a result of the m = 2 component of the tidal field 
(see eq. @). 

In several X-ray binaries, most notably Her X-l, there is evidence for a tilted, precessing disk (see Wijers & Pringlc 
1999 and references therein). The tilting mechanism we have described is very unlikely to operate in such disks, which are 
expected to be tidally truncated and to have H/r <c 0.4. Therefore, it appears that tidal torques are not responsible for 
the tilting of disks in X-ray binaries (cf. Larwood 1998). Possible mechanisms for tilting these disks include wind torques 
(Schandl & Meyer 1994) and radiation torques (Wijers & Pringle 1999). 

Another possible application of this work is to nearly Keplerian disks that surround black holes in active galactic nuclei. 
If the disk is subject to a bar potential from the galaxy and the disk radius is sufficiently smaller than the corotation 
radius of the bar, then the disk will be subject to this tilt instability. 

The results in this paper have implications to protostellar disks perturbed by inclined planets. A secular resonance 
occurs where the precession frequency of a planet matches the local precession frequency of an orbiting particle. The 
resonant radius changes as the nebula disperses and the resonance sweeps across a major portion of the solar nebula 
(Ward 1981). However, the current results suggest that the effects of such resonances on the gaseous nebula are mild and 
are distributed over the disk. Further analysis can be carried out through the methods described in this paper. 

We thank Jim Pringle for encouraging this investigation and for providing useful discussions. We acknowledge support 
from NASA grant NAG5-4310 and from the STScI visitor program. GIO was supported by the European Commission 
through the TMR network 'Accretion on to Black Holes, Compact Stars and Protostars' (contract number ERBFMRX- 
CT98-0195). 
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Fig. 1. — Tilt instability cycle induced by the binary tidal field. We label disturbances by integer pairs (m, /) for azimuthal and temporal 
dependences in the inertial frame of the form exp[i(— m<f) + IQ^t)], for azimuthal angle (p. The tilt and tide interact to produce a bending 
wave. The wave and tide interact to produce a stress that amplifies the tilt. 



Lubow & Ogilvie 



15 




0.0008 
0.0006 

a 0.0004 

? 0.0002 

0.0000 

-0.0002 

0.0 




0.1 0.2 



0.3 



0.004 
0.003 

a 0.002 

? 0.001 




0.3 



Fig. 2. — Growth rate of the modified rigid-tilt mode, in units of the binary frequency, plotted against the outer radius of the disk, in units 
of the binary radius. The four panels are for a = 0.1 (top left), a = 0.01 (top right), a = 0.001 (bottom left), and a = (bottom right); 
otherwise, parameters have their reference values. Note that the vertical scale is different in each case. 
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Fig. 3. — Expanded view of Fig. 2, showing the primary resonance in the cases a = (solid line) and a = 0.001 (dotted line). 
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Fig. 4. — Precession frequency of the modified rigid-tilt mode, in units of the binary frequency, plotted against the outer radius of the disk, 
in units of the binary radius. Solid line: reference model (a = 0.01). Dashed line: analytic approximation from Bate et al. (2000). 
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Fig. 5. — Decay rate of the modified rigid-tilt mode, in units of the binary frequency, divided by a and plotted against the outer radius 
of the disk, in units of the binary radius. Solid line: reference model (a = 0.01). Dotted line: reference model, but with only the m = 
component of the tidal potential. Dashed line: analytic approximation from Bate et al. (2000). 
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(al) (a2) (a3) (a4) 




Fig. 6. — Shape of the disk while executing the modified rigid-tilt mode. The disk has r2/Vb = 0.118, in the middle of the primary resonance. 
Panels (a), (b), and (c) correspond to the cases a = 0, a = 0.001, and a = 0.01, respectively, with growth rates — uii/Q^, = 0.003596, 0.000655, 
and 0.000065. In each case, view (1) is an :rz-scction, looking along the j/-axis, at phase of the cycle seen in the binary frame. View (2) is 
a j/z-section, looking along the negative rr-axis, as if from the companion star. Panels (3) and (4) are the same as (1) and (2), but at phase 
7r/2 in the cycle. The eigenfunctions have been renormalized for ease of comparison. 




Fig. 7. — As for Fig. 6, but for a disk with r2/Vb = 0.3. Here a = 0.01 is fixed and panels (a), (b), and (c) correspond to the cases e = 0.1, 
e = 0.05, and e = 0.03, respectively, with damping rates uj { /n h = 0.000258, 0.001205, and 0.003601. 
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APPENDIX 

DERIVATION OF THE REDUCED EQUATIONS FOR LINEAR BENDING WAVES 
Let the small parameter e be a characteristic value of the angular semi-thickness H/r of the disk. Then set 

K 2 = n 2 [l + ef K (r)}, (Al) 
tf z = tf[l + ef z (r)}, (A2) 
a = ef a (r), (A3) 

(A4) 

where the functions / are O(l), which includes the possibility of their being arbitrarily small. The disk is assumed to 
satisfy the Navier-Stokes equation with a (dynamic) shear viscosity given by 

<«) 

To describe the internal structure of the disk, adopt units in which the radius of the disk and the orbital frequency are 
0(1). Introduce the stretched vertical coordinate 

C = ~ e , (A6) 

which is O(l) inside the disk. We then find, for the unperturbed disk, 

U = 0(e 3 ), (A7) 

v = rn(r)+e 2 rn 2 (r,()+0(e 3 ), (A8) 

w = 0(e 4 ), (A9) 

p = e s [p (r,C) + e Pl (r,0 + O(e 2 )], (AlO) 

P = e s+2 [ Po (r, C) + e Pl (r, () + 0(e 2 )] , (All) 

A1 = e s + 3 Mr,C) + 0(e)], (A12) 

where s is an arbitrary positive parameter. Any viscous evolution of the disk occurs on a long time-scale 0(e~ 3 ) and is 
consistently neglected. The vertical component of the equation of motion implies, at O(e), 

= --|^ 2 C, (A13) 



and, at 0(e 2 ), 

while the radial component at 0(e 2 ) gives 



1 dpi pi dpg 
Po d( pi d( 



+ ^^"/^ 2 C (AW) 



2H^ 2 = -1^-4V- (A15) 
po or dr 



Consider linear bending waves with azimuthal wavenumber m = 1 in which the Eulerian perturbation of u, say, is 

Re [u'(r, z, t) e - **] . (A16) 

It is known that these waves travel radially at a speed comparable to the sound speed (Papaloizou & Lin 1995). Therefore 
the characteristic time-scale for the evolution of the warped shape is the radial sound crossing time r/c s « e _1 ri _1 , 
implying that the perturbations evolve on a time-scale 0(e~ 1 ) that is long compared to the orbital time-scale [0(1)] but 
much shorter than the viscous time-scale [0(e -3 )]. This is captured by a slow time coordinate 



et. (A17) 



For the perturbations, introduce the sea" 



mgs 

= eu[(r, C, T) + e 2 u' 2 (r, £ T) + 0(e 3 ), (A18) 

= ev[(r, C, T) + e 2 ^(r, C, T) + 0(e 3 ), (A19) 

= ew[(r, C, T) + e 2 w' 2 (r, (, T) + 0(e 3 ), (A20) 

= e s [p[(r, C, T) + ep' 2 (r, (, T) + 0(e 2 )] , (A21) 

p' = e s + 2 [p'^r, C, T) + ep' 2 (r, (, T) + 0(e 2 )] . (A22) 
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The overall amplitude of the perturbations is of course arbitrary since this is a linear analysis. |^] 
The perturbed equations for w, p, and p at leading order are 



■o / , 1 d Pi Pi d Po n , AO o A 
illw 1 + — 2 -— = 0, (A23) 

Po d( pi d( 



where 7 is the adiabatic exponent. These may be combined to give 



- i np[ + w> 1 ^+p o ^ = 0, (A24) 
iOpi + w[ % + 7 Po ^ = 0, ( A25) 



The general solution, regular at the disk surface, is 



If^hO. (A26) 



1 



: irQW, (A27) 

p[=r^W, (A28) 

p[=r^W, (A29) 

where W(r, T) is a dimensionless complex function to be determined. These perturbations correspond to applying a rigid 

tilt to each annulus of the disk. The tilt varies with radius and time according to the function W{r, T), which is related 
to the unit tilt vector £ through 

W = t x + i£ y . (A30) 
The perturbed equations for u and v at leading order are 

- iQu[ - 2Qv[ = 0, (A31) 

-iQv[ + ifiu; = 0. (A32) 



The general solution is 



u[ = U, (A33) 
v[ = -\iU, (A34) 



where U[r, C, T) is a complex function to be determined. 

The perturbed equations for w, p, and p at the next order are 



1 dp' 2 _ p'i dpa 
Po dC pi dC, 



inw ' 2 + ±^.-q^=F w , (A35) 



where 



-m P > 2 +w^ + p ^ = F p , (A36) 
-m P > 2 +w' 2 ^ + 7 Po^ = F p , (A37) 

Fw = -zrtt^ - ^rtfW - ^f z rtf(W, (A38) 
oT p Q p Q dC 

dpo dW . n dpi Id p U 

ac ~d¥ ~ irQ ~dc w ~ rYr {porU) + 

5 The scaling adopted here is, however, of some significance, since it corresponds to a (differential) tilt angle comparable to the angular 
thickness of the disk. This is appropriate for a warp of observational consequence. We note that nonlinear effects may be significant for warps 
of this amplitude. Unfortunately, the Eulerian perturbation method used here tends to overestimate the degree of nonlinearity. For example, 
although the fractional Eulerian density perturbation is of order unity, the dominant perturbation is (locally) a rigid tilt and the Lagrangian 
density perturbation is in fact of higher order in t. Nonlinear effects can occur, however, owing to the fact that the horizontal motions are 
comparable to the sound speed. In particular, these motions can be damped by a parametric instability (Gammie, Goodman, & Ogilvie 
2000; Bate et al. 2000). The Eulerian method remains the most convenient way of obtaining the equations if one is satisfied with a formal 
linearization. Since we are considering a stability problem in this paper, this method is sufficient for our purposes. 



F P = ~r-^^F- irQ^W - ~(p rU) + ^, (A39) 
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and F p will not be required. Now the linear operator defined by the left-hand sides of equations ( A35 )-( A37 ) is singular 
owing to the existence of the tilt mode identified above. The corresponding solvability condition is 



J (i Po F w + FpClQ d( = 0, 



where the integral is over the entire vertical extent of the disk. This evaluates to 



dW 



1 d 



where 



^,r<> i 2^ +if z nWj - ^ / p r 2 W(<!^ = 0. 
S = po d( 



is the surface density. 

The perturbed equations for u and v at the next order are 



iQu' 2 — 2Qv' 2 = F Ul 
-iClv' 2 + 4\nu' 2 = F v 



where 



dU 19/ dp. 



dT po dr V dC, 
dU 



r dp dp . 



1 d 



pl d( dr p d( \ d( 



dU 



dn 



i dpo. 



d ( dU 



dT d( po d( 2p dC \ d( 

Again the linear operator is singular, with the solvability condition 



This evaluates to 



F u + 2iF v = 
dU 



9 9U -if OTT+ 2 f<* 9 

or Ik u + Po n d( y u d( 



dW 



po— ) +m z (— = o 



By inspection, the solution is of the form U oc £, with 

d 



2— + if K n + 2f a n - = m 2 — . 



dT 



u 



,dW 



dr 



If we now define 



where 



G = 



l r 2 n fU 

c 



Zo = / Pot 2 d( 

is the second vertical moment of the density, we obtain the coupled equations 



Finally, if we step back from the asymptotic analysis and present the equations in physical terms, we obtain 



Er 2 n 

dG 
dt ' 



dw (nl -n 2 \ in 



dt 



n 2 



-w 



n 2 \ in 



n 2 



ldG 

r dr ' 

Tr 3 n 3 dW 



where now 



—G + aQG = 4 



S = / pdz, 



pz dz. 



(A40) 

(A41) 
(A42) 



(A43) 
(A44) 



(A45) 
(A46) 

(A47) 
(A48) 

(A49) 

(A50) 
(A51) 

(A52) 
(A53) 

(A54) 
(A55) 

(A56) 
(A57) 



